function [iteration_number,convergent_fraction,convergent_fraction_b, ...
    divergent_fraction,divergent_fraction_b] ...
    = plot_histograms(n_iterations,n_iterations_b,n_realizations)

% Function to compute fractions of convergence and divergence of the
% Newton-Raphson scheme, and to generate histograms (convergence with
% respect to iteration number)

% Compute gross fractions of convergence/divergence
n_convergent = n_iterations(~isnan(n_iterations));
n_convergent_b = n_iterations_b(~isnan(n_iterations_b));
divergent_fraction = 1 - numel(n_convergent)/n_realizations;
divergent_fraction_b = 1 - numel(n_convergent_b)/n_realizations;


% Compute histogram (fraction which converges at each iteration value)
convergent_min = min(n_convergent);
convergent_max = max(n_convergent);
convergent_b_min = min(n_convergent_b);
convergent_b_max = max(n_convergent_b);
if isempty(convergent_b_min)
    min_hist = convergent_min;
else
    min_hist = min(convergent_min,convergent_b_min);
end
if isempty(convergent_b_max)
    max_hist = convergent_max;
else
    max_hist = max(convergent_max,convergent_b_max);
end
for ii = min_hist: max_hist
    iteration_number(ii-min_hist+1) = ii;
    convergent_fraction(ii-min_hist+1) = sum(n_convergent==ii)/n_realizations;
    convergent_fraction_b(ii-min_hist+1) = sum(n_convergent_b==ii)/n_realizations;
end
% plot histograms
figure(100)
hold on
plot(iteration_number,convergent_fraction,'ko','MarkerSize',7)
plot(iteration_number,convergent_fraction_b,'kx','MarkerSize',9)
xlabel('Number of iterations')
ylabel('Convergent fraction')
xticks(iteration_number);
ax = gca; ax.FontSize = 16;
% figure out good limits to put on plot
xlim([iteration_number(1)-2, iteration_number(end)+5])
y_max_convergent = max(max(convergent_fraction,convergent_fraction_b));
y_max_divergent = max(divergent_fraction,divergent_fraction_b);
y_max = max(y_max_convergent,y_max_divergent);
ylim([-0.05, y_max+0.05])
% plot divergent proportion to the right
divergent_location = iteration_number(end)+4;
xline(divergent_location-1,'-',{'Divergent Fraction'},'LineWidth',2);
plot(divergent_location,divergent_fraction,'ko','MarkerSize',7)
plot(divergent_location,divergent_fraction_b,'kx','MarkerSize',9)


